In kmcd39/divM: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(echo = TRUE)

rm(list=ls())
library(sf)
library(dplyr)
library(purrr)
library(mapview)
library(lwgeom)
#library(divM)
devtools::load_all()

This markdown walks through the generation of different versions of the polygon dividedness measure with the NHPN highway data. First it walks through different ways to setup/subset the NHPN highway data, and then shows how these affect the number of polygon subdivisions generated by the measure.

Load region + division boundaries to begin:

# load & setup NHPN hwy data
# https://catalog.data.gov/dataset/national-highway-planning-network-nhpn
shp.dir <- "~/R/shapefiles/" %>% 
hwys <- st_read(paste0(shp.dir, "National_Highway_Planning_Network-shp/National_Highway_Planning_Network.shp"))

# get czs from my data storage pkg
czs <- divDat::czs

# make uniform metered crs -----------------------------------------------
czs <- czs %>% conic.transform()
hwys <- hwys %>% conic.transform()

Subsetting by hwy types

There are many different ways to subset hwys within a region for generating the measures:

• By route type: "Interstate," "State route," "non-signed other arterial," etc.. (SIGNT1 in data)
• By "Functional classification": --- can get to "limited-access" hwys this way. (FCLASS in data)
• By spatial relationship: Interstates and all hwys (of X types) that intersect interstates

Illustration case for philly:

# get sample area & trim hwys
(phl <- czs[grepl("Philadelphia",czs$cz_name),] %>% divM::region.reorg("cz")) 
hwys <- hwys %>% st_crop(phl)

# use FCLASS to create a hwy subset
lac <- hwys %>% filter(FCLASS %in% divM::lac_codes |
                          SIGNT1 == "I") # all limited-access

# spatial group hwy data
hwys <- hwys %>% denode.lines(group.cols = c("SIGNT1", "SIGN1"))
lac <- lac %>% denode.lines(group.cols = c("SIGNT1", "SIGN1"))

# Get hwy subsets for philadelphia: ----
# all hwys
all <- subset.polys.divs(region = phl,
                               div.sf = hwys)
# all interstates
int <- subset.polys.divs(phl, hwys,
                               div.identifier.column = "SIGNT1",
                               always.include = "I", 
                               include.intersecting = F) 
 # all interstates and hwys that intersect interstates
tint <- subset.polys.divs(phl,
                                hwys,
                                div.identifier.column = "SIGNT1",
                                always.include = "I",
                                include.intersecting = T)


lac <- subset.polys.divs(phl,
                                lac)

# visualize differences
plot(all['SIGN1'], main = "all philly hwys", lwd=2)
plot(int['SIGN1'], main = "all philly interstates", lwd=2)
plot(tint['SIGN1'], main = "all philly interstates+touching",lwd=2)
plot(lac['SIGN1'], main = "all philly LAC",lwd=2)

"Filling gaps" / The linebreak problem

Sometimes a hwy has a minuscule break. Often the hwy classification changes for an interchange or a short stretch. For example, I676 in Camden.

I run a helper fcn to fill these before generating polygons. The logic of the function is: If a line segment has an endpoint within x threshold of another segment, and they have the same hwy identifiers, I connect the segments.

Default threshold for filling gap is 200 meters. The argument threshold can be passed to any wrapper function that fills gaps to change this.

Below are the 3 interstates in Philadelphia CZ that have these gaps.

# filling gaps in Philly interstates:
viz.fix <- Fix.all.hwys(int, return.gap.map = T, verbose = T)

viz.fix$I676
viz.fix$I76
viz.fix$I95

To continue creating the polygon measures, create versions of the above hwy subsets that have their gaps filled. Suffix names with '.c'

int.c <- suppressWarnings(Fix.all.hwys(int, verbose = F))
tint.c <- suppressWarnings(Fix.all.hwys(tint, verbose = F))
lac.c <- suppressWarnings(Fix.all.hwys(lac, verbose = F))

Creating polygon measures

From here, creating the polygon measure is pretty straightforward. Fcn polygonal.div does the work; it takes linear division data and a polygon representing study area.

Let's see how it looks based on the different hwy subsets and whether or not the gaps are filled. I use CZ as boundary.

The other core parameter for generating these measures is the area floor--- the minimum size of each polygon subdivisions. If a polygon is formed that is below this area, it is filtered out and not counted towards the total. By default this floor is 1/2 square km.

# Summary table:
Summary.table <- 
  tibble( gaps_filled = c(FALSE, TRUE)
           ,interstates_only = c(polygonal.div(phl, int)$n.polys,
                                 polygonal.div(phl, int.c)$n.polys)
           ,interstates_intersecting = c(polygonal.div(phl, tint)$n.polys,
                                         polygonal.div(phl, tint.c)$n.polys)
           ,limited_access = c(polygonal.div(phl, lac)$n.polys,
                                 polygonal.div(phl, lac.c)$n.polys)
           )

knitr::kable(Summary.table)

# to view the polygons as map, ask the same function to return.sf to map
polygonal.div(phl, lac.c, return.sf = T) %>%
  select("id") %>%
  plot(main = "LAC polygons in philly")

Wrapping Up

To end on a pun, a single function, Polys.wrapper, can replicate this whole workflow: of subsetting the hwys, filling the gaps, and generating the polygons. It will pass on any arguments to its wrapped fcns. See full options with ?Polys.wrapper.

suppressWarnings( Polys.wrapper(phl, lac,
              fill.gaps = T, verbose = F) )

kmcd39/divM documentation built on Oct. 21, 2023, 11:28 p.m.